Matrix approach to discrete fractional calculus II: 
partial fractional differential equations 



Igor Podlubny^'*, Aleksei Chechkin'', Tomas Skovranek^, YangQuan Chen^, 

Bias M. Vinagre Jara*^ 

"■BERG Faculty. Technical University of Kosice, Slovak Republic 
^Institute for Theoretical Physics NSC KIPT, Kharkov, Ukraine 
'^Department of Electrical and Computer Engineering, Utah State University, USA 
'''Industrial Engineering School, University of Extremadura, Badajoz, Spain 



Abstract 

A new method that enables easy and convenient discretization of partial 
differential equations with derivatives of arbitrary real order (so-called frac- 
tional derivatives) and delays is presented and illustrated on numerical solu- 
tion of various types of fractional diffusion equation. The suggested method 
is the development of Podlubny's matrix approach (Fractional Calculus and 
Applied Analysis, vol. 3, no. 4, 2000, 359-386). Four examples of nu- 
merical solution of fractional diffusion equation with various combinations 
of time/space fractional derivatives (integer/integer, fractional/integer, in- 
teger/fractional, and fractional/fractional) with respect to time and to the 
spatial variable are provided in order to illustrate how simple and general is 
the suggested approach. The fifth example illustrates that the method can 
be equally simply used for fractional differential equations with delays. A 
set of MATLAB routines for the implementation of the method as well as 
sample code used to solve the examples have been developed. 
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1. Introduction 

Recently, kinetic equations of the diffusion, diffusion-advection, and 
Fokker-Planck type with partial fractional derivatives were recognized as 
a useful approach for the description of transport dynamics in complex sys- 
tems whose temporal evolution deviates from the standard laws, that is, 
from exponential Debye or Gaussian laws, and from fast decaying correla- 
tions. Examples include systems exhibiting Hamiltonian chaos, disordered 
medium, plasma and fluid turbulence, underground water pollution, dynam- 
ics of protein molecules, motions under the influence of optical tweezers, 
reactions in complex systems, and more (see reviews on fractional kinet- 
ics [HI HH Uni EH US], the recent multi-author book [28], and references 
therein). These fractional equations are derived asymptotically from basic 
random walk models, the generalized master and Langevin equations. The 
advantage of the fractional models lies in the straightforward way of including 
external force terms and of calculating boundary value problems. Also, the 
consideration of transport in the phase space spanned by both position and 
velocity coordinates is possible within the fractional approach. However, be- 
cause of complicated integro-differential structure of fractional kinetic equa- 
tions the analytical solutions are presently known only in a very few relatively 
simple cases. Therefore, the development of numerical methods is of current 
importance. 

Let us recall briefly how the kinetic equations with integer partial deriva- 
tives can be "fractionalized" . There are two generic types of fractionaliza- 
tion, which can be explained by taking as an example the parabolic diffusion 
equation for the particles density u{x, t) in a one-dimensional space, 

^ = X^, (t>0, a<a;<6) (1) 

where constant x is diffusion coefficient. The first type of fractionalization 
leads to a time fractional diffusion equation by means of replacing the first 
order time derivative by afractional derivative of order a less than 1, 

^D^u = Xg^, {t>0, a<x<b) (2) 
Here, qD" is the Caputo fractional derivative [2], which is defined as 
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^D^Jix) = ^^_ ^ / ^7 , (m-l</i<m) (3) 



{m- fi) J (x - ^ 

a 

Taking a = 1 in ([2]) gives the classical diffusion equation ([T]). 

Other two forms of a time fractional diffusion equation that appears in 
the literature use the Riemann-Liouville fractional derivative instead of the 
Caputo one [H]. Although recently, in addition to a geometric and physical 
interpretation of fractional integration and fractional differentiation [52j, a 
physical interpretation for the initial conditions in terms of the Riemann- 
Liouville fractional derivatives of the unknown function has been suggested 
[27] . the use of Caputo derivative in physical problems is perhaps more conve- 
nient since it allows using initial conditions expressed in terms of values of the 
unknown function and its integer-order derivatives [HO]. However, all three 
forms of "time-fractionalization" are equivalent if zero initial conditions are 
posed. In what follows we use the form with the Caputo derivative, equation 
(|2]), since some of the illustrating examples use non-zero initial conditions. 

In the second type of fractionalization, the second order spatial derivative 
is replaced by the fractional derivative of the order (3 between 1 and 2, thus 
leading to spatial fractional diffusion equation, 

du d^u , , , . , 

^=X^, (t>0, a<.<h) (4) 

where ld\x\^ (we adopt here the notation introduced in [M]) is a partial 
(with respect to spatial variable) symmetric Riesz derivative, which is defined 
as a half-sum of the left- and right-sided Riemann-Liouville derivatives [501 
EI]: 



d\x\ 



Dl<^{x) = -[aBim^ ^Dl<^{x)), (5) 



where the left- and right-sided Riemann-Liouville derivatives are defined by 
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xD'rcbix) = — r — — / — ^ — , m — 1 < u < m , 7 

''^^ ' r(m-/i) V dxj J ^ p_ 

X 

For (3 = 2 the equation Q becomes the classical diffusion equation ([T]). 

Other forms of asymmetric space fractional generalizations use the left- 
side Riemann - Liouville derivative instead of the symmetric Riesz deriva- 
tive [12l HO], or asymmetric derivative with different asymmetry factors at 
the left- and right-side derivatives [THl [TTl |32]- In terms of random walk 
schemes, the symmetric derivative corresponds to a symmetric jump prob- 
ability distribution of a diffusing particle, whereas any asymmetry in space 
derivative accounts for inherent force-free preferable direction of jumps which 
may occur, e.g., in heterogeneous porous media or magnetically confined fu- 
sion plasmas. In our paper we restrict ourselves to symmetric case, equation 

Of course, there are different generalizations of time and space fractional 
diffusion equations, including: multidimensional fractional diffusion and ki- 
netic equations [5l|T7l|38], both time and space fractional generalizations [36] , 
different regular forces in space and time fractional Fokker-Planck equa- 
tions in [m E?! El], variable transport coefficients [B3], equations 
with fractional derivatives of distributed and variable orders [TJ El El EZj 
etc. The realm of fractional kinetics is growing, and therefore it is desirable 
to have at hand a method for numerical solution which would be relatively 
simple and at the same time general enough to deal effectively with differ- 
ent forms of fractional kinetic equations. However, while different numerical 
tools for ordinary fractional equations exist and a basic framework of their 
numerical solution is already established, relatively few numerical methods 
exist to solve fractional equations with partial derivatives, and the develop- 
ment of effective numerical schemes is now on the agenda. We recall briefly 
the different approaches used in the literature. 

The numerical methods differ essentially in the way in which normal and 
fractional derivatives are discretized. In [35] to solve diffusion-reaction equa- 
tion with the left Riemann-Liouville derivative between 1 and 2, the L2 dis- 
cretization method was used taken from [17], together with its modification, 
L2C (both L2 and L2C methods are based on numerical approximation of 
a fractional integral that appears in the definition of the Riemann-Liouville 
fractional derivative). It was shown that the former is the most accurate 
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for orders larger than 1.5, whereas the latter is the most accurate for orders 
less than 1.5. For the first order time derivative, the explicit forward Euler 
formula and semi-implicit scheme were used. 

Langlands and Henry |29] used LI scheme from ^7] to discretize the 
Riemann-Liouville fractional time derivative of order between 1 and 2. 

Yuste [SZ] considered a Griinwald-Letnikov approximation for the Riemann- 
Liouville time derivative and used a weighted average for the second-order 
space derivative. 

Scherer et al. [55] introduced very recently a modification of the Griinwald- 
Letnikov approximation for the case of the Caputo derivative of a function 
which is not zero in the starting point of the considered time interval, and 
applied that approximation for the numerical solution of fractional diffusion 
equations with the Caputo time derivative and non-zero initial conditions. 

To solve the one-dimensional space fractional advection-dispersion equa- 
tion with left-side Riemann-Liouville derivative and variable coefficients the 
shifted Griinwald-Letnikov approximation was proposed by Meerschaert and 
Tadjeran [ID]. For two-sided space-fractional partial differential equations the 
shifted Griinwald-Letnikov formula was proposed and discussed in [H]. The 
fractional Crank-Nicholson method based on the shifted formula was elabo- 
rated, giving temporally and spatially second-order numerical estimates [61] . 
The generalizations of the shifted formula and of the fractional Crank- Nicholson 
method in the two-dimensional case were discussed in [32] and [SO], respec- 
tively. 

Another method to solve the space-fractional Fokker-Planck equation 
with constant coefficient on the fractional derivative term was pursued by 
Liu et al. [32] • They transform the partial differential equation into a sys- 
tem of ordinary differential equations, which is solved by a method of lines. 

Ervin and Roop [151 [IB] presented a theoretical framework for the Galerkin 
finite element approximation to the steady state fractional advection-diffusion 
equation, and extended this approach to multidimensional partial differential 
equations with constant coefficients on the fractional derivative terms. 

Valko and Abate [^2] solved the time-fractional diffusion equation on a 
semi-infinite domain by numerical inversion of the two-dimensional Laplace 
transform. To solve the time-fractional diffusion equation in a bounded do- 
main, Lin and Xu [31] proposed the method based on a finite difference 
scheme in time and Legendre spectral method in space. 

Liang and Chen [30] used a combination of symbolic computations and 
numerical inversion of the Laplace transform for solving a time-fractional 
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diffusion-wave equation with the time derivative of order between 1 and 2. 

We also mention that in order to approximate shifted Caputo time deriva- 
tive appearing in hydrodynamic equations for heterogeneous porous media 
the modification of Yuan and Agrawal's method |66] was used to transform a 
fractional derivative into an infinite integral over auxiliary internal variables 

Another approach for the solution of fractional kinetic equations employs 
the methods of Monte Carlo type (random walk based methods). A set of 
random walk schemes applied to fractional diffusion equations based on the 
Griinwald-Letnikov approximation was developed in the papers by Goren- 
fio, Mainardi and co-workers. They were applied to solve (i) symmetric 
space- fractional diffusion equation [20l |22] ; (ii) asymmetric space-fractional 
diffusion equation in the Levy-Feller form [21]; (iii) time-fractional diffu- 
sion equation with Caputo derivative [25]; (iv) time-space fractional diffu- 
sion equation [2H |23]. Chechkin et al. [9] generalized the approach on a 
double-order time fractional diffusion equation. Gorenflo and Abdel-Rehim 
[19] proposed discrete approximations to time-fractional diffusion process 
with non-homogeneous drift towards the origin by generalization of Ehren- 
fest's urn model. The Levy-Feller diffusion-advection process with a con- 
stant drift was approximated by random walk and finite difference method 
by Liu et al. [23] • The random walk particle tracking approach to solve 
one-dimensional space-fractional diffusion-advection equation with space de- 
pendent coefficients was employed by Meerschaert and co-authors [65]. The 
method based on numerical solution of a coupled stochastic differential equa- 
tions driven by Levy symmetric stable processes was proposed in [58] to solve 
a non-linear evolution problem involving the fractional Laplacian operator. 

All aforementioned works indicate that numerical solution of partial frac- 
tional differential equations plays an important and increasing role in the 
applications of the methods and models of non-integer order. 

In the present paper we propose a general approach to the numerical 
solution of partial fractional differential equations, which is based on the 
matrix form representation of discretized fractional operators introduced in 
|51j . This approach unifies the numerical differentiation of arbitrary (includ- 
ing integer) order and the n-fold integration, using the so-called triangular 
matrices. Applied to numerical solution of differential equations, it also uni- 
fies the solution of integer- and fractional-order partial differential equations. 
The suggested approach leads to significant simplification of the numerical 
solution of partial differential equations, and it is general enough to deal with 
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different types of partial fractional differential equations, even with delays. 

2. The idea of the suggested method 

The method that we suggest is based on triangular strip matrix approach 
[5T] to discretization of operators of differentiation and integration of arbi- 
trary real order. 

In contrast with generally used numerical methods, where the solution 
is obtained step-by-step by moving from the previous time layer to the next 
one, let us consider the whole time interval of interest at once. This allows 
us to create a net of discretization nodes. In the simplest case of one spatial 
dimension this step gives a 2D net of nodes. An example of such discretization 
is shown in Fig. [l] The values of the unknown function in inner nodes (shaded 
area in Fig. [T| are to be found. The values at the boundaries are known: 
they are used later in constructing the system of algebraic equations. 

The system of algebraic equations is obtained by approximating the equa- 
tion in all inner nodes simultaneously (this gives the left-hand side of the 
resulting system of algebraic equations) and then utilizing the initial and 
boundary conditions (the values of which appear in the right-hand side of 
the resulting system). 

The discretization nodes in Fig. [T] are numbered from right to left in each 
time level, and the time levels are numbered from bottom to top. We use 
such numbering in this article for the clarity of presentation of our approach, 
although standard numberings work equally well. 

In the following sections we recall the basic tools that are necessary for 
the method: the triangular strip matrices, the Kronecker product, the elimi- 
nators, and the shifters. Then we show how they are used for approximating 
partial derivatives of arbitrary real order and the equation, and how the 
resulting system of algebraic equations appears. 

3. Triangular strip matrices 

In this paper we use matrices of a specific structure, which are called 
triangular strip matrices [SH EH], and which have been also mentioned in 
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Figure 1: Nodes and their right-to- left, and bottom-to-top number 



[H HE]. We will need lower triangular strip matrices, 
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A lower (upper) triangular strip matrix is completely described by its first 
column (row). Therefore, if we define the truncation operation, truncAr(-), 
which truncates (in a general case) the power series g{z), 



fc=0 



(10) 



to the polynomial Qn^z) 



N 



truncAT {g{z)) '= 



iOkZ 



Qn{z), 



(11) 



k=0 



then we can consider the function g{z) as a generating series for the set of 
lower (or upper) triangular matrices (or Un), N = 1, 2, . . . 

It was shown in |51] that operations with triangular strip matrices, such 
as addition, subtraction, multiplication, and inversion, can be expressed in 



the form of operations with their generating series (10). 



Among properties of triangular strip matrices it should be noticed that 
if matrices C and D are both lower (upper) triangular strip matrices, then 
they commute: 

CD = DC. (12) 



4. Kronecker matrix product 

The Kronecker product A^Boi the nxm matrix A and the pxq matrix 



B 



A 



ail 0-12 • • • Cllm 
021 0-22 ■ ■ ■ 0,2m 

Onl On2 ■ ■ ■ OjYim 



B 



bii bi2 ... big 

621 &22 • • • b 



2q 



bpi bp2 ... b 



pq 



is the np x mq matrix having the following block structure: 



A®B 



aii-B ai2-B . . . aimB 
021B a22B . . . a2mB 

ttnlB a„2-B . . . CLnmB 



(13) 



(14) 
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For example, if 
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B 



1 2 

4 5 
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(15) 



(16) 



Among many known interesting properties of the Kronecker product we 
would like to recall those that are important for the subsequent sections. 
Namely [63j, 

• ii A and B are band matrices, then A® B also a band matrix, 

• ii A and B are lower (upper) triangular, then A®B is also lower (upper) 
triangular. 

We will also need two specific Kronecker products, namely the products 
En® A and A ® Em-, where is an n x n identity matrix. For example, if 
A is a 2 X 3 matrix 
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This illustrates that left multiplication of Anxm by En creates an n x n 
block diagonal matrix by repeating the matrix A on the diagonal, and that 
right multiplication of Anxm by Em creates a sparse matrix made of n x m 
diagonal blocks. 



5. Eliminators 



The suggested method requires also the use of a certain type of matrices 
called eliminators , which are obtained from the N x N unit matrix E by 
keeping only some of its rows and omitting all other rows: 5*1 is obtained by 
omitting only the first row of E; S2 is obtained by omitting only the second 
row; Si^2 is obtained by omitting only the first and the second row of E; 
and, in general, Sri,r2,...,rk is obtained by omitting the rows with the numbers 
ri, r2, . . . , rfc. In case of infinite matrices, similar matrices appeared in [TU] . 

If A is a square N x N matrix, then the product Sr-i^^r2,...,rk^ contains 
only rows of A with the numbers different from ri, r2, . . . , r^. Similarly, the 
product AS"^ j.^^ contains only columns of A with the numbers different 
from ri,r2, . . . ,rfc. 

The following simple example illustrates the main property of eliminators: 
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0,22 O23 
O32 '^33 



O2I O22 ^23 
O3I 0,32 ^33 



6. Shifters 

For some types of approximation of differential operators (for example, 
one of the approximations of the symmetric Riesz derivative below in this 
article) and especially for numerical solution of differential equations of ar- 
bitrary order (integer or fractional) with delays, it is convenient to introduce 
another special kind of matrices - shifters -, which will represent discrete 
shifts, like, for example, delays. 

Shifters (although without using this term) were used in [51] for a simple 
generation of triangular strip matrices. There are shifters of two kinds: (A^ + 
1) X (A^ + 1) matrices E^^, p = 1, . . . N, with ones on p-th diagonal above the 
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main diagonal and zeroes elsewhere, and matrices Ej^^, p = 1, . . . N, with 
ones on p-th diagonal below the main diagonal and zeroes elsewhere. We 
also denote = En the unit matrix. 

The shift of all the coefficients in the triangular strip matrix f/^r in the 
south-west (bottom-left) direction can be easily written if we start with Un+i 
and then use shifters and eliminators: 

-iUn = Si Ej^^^^^ Un+1 -E'iv+i,! 'S'^+i (20) 

Similarly, the shift of all the coefficients in the triangular strip matrix Un 
in the north-east (top-right) direction can be easily obtained as: 

+iUn = Sn+1 E'}^_^_^^^Un+i E^_^_^^^ Si (21) 



7. Discretization of ordinary fractional derivatives 

It follows from [5T], that the left-sided Riemann-Liouville or Caputo frac- 
tional derivative v^"\t) = QD^v{t) can be approximated in all nodes of the 
equidistant discretization net t = jr {j = 0,1, ... ,n) simultaneously with 
the help of the upper triangular strip matrix i^n""* as [_ : 



where 



-'n-l 



V 



1 T 



UJ, 

















UJ 







1 



UJ 



OJi 







Vn Vn-1 



1 



UJ 








UJ. 
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(a) 
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(a) 
J 



(22) 



(23) 



^ In this article due to the use of the descending numbering of discretization nodes 
the roles of the matrices i?!"' (originally for backward fractional differences) and F^"^ 
(originally for forward fractional differences) are swapped in comparison with |51j . where 
these matrices were introduced for the first time. However, we would like to preserve the 
notation bIi'^ for the case of the backward fractional differences approximation and -F^") 
for the case of the forward fractional differences approximation. 
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(a) 



-11 



3 = 0, 1, 



(24) 



Similarly, the right-sided Riemann-Liouville or Caputo fractional deriva- 
tive = tDb^it) can be approximated in all nodes of the equidistant 
discretization net t = jr (j = 0, 1, . . . , n) simultaneously with the help of the 
lower triangular strip matrix Fn 
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(25) 



UJ 
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UJ. 



n — 1 



UJ, 



UJ 



UJ, 



(a) 



(26) 



The symmetric Riesz derivative of order f3 can be approximated based 
on its definition (|5]) as a combination of the approximations (22) and (25) 
for the left- and right-sided Riemann-Liouville derivatives, or using the cen- 
tred fractional differences approximation of the symmetric Riesz derivative 
suggested recently by Ortigueira [Ml US]. The general formula is the same: 



-^m—l 



V 







^m— 1 



Vi Vq 



(27) 



In the first case, the approximation for the left-sided Caputo derivative 
is taken one step ahead, and the approximation for the right-sided Caputo 
derivatve is taken one step back. This leads to the matrix 



(28) 



In the second case (Ortigueira's definition |18]) we have the following 
symmetric matix: 
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(29) 



W) 



;-l)^r(/3 + l) cos(/37r/2) 



r(/5/2-fc + l)r(/?/2 + A; + l)' ^-O'l'---'"^ (30) 

Both these approximations of symmetric Riesz derivatives give practi- 
cally the same numerical results and in case of numerical solution of partial 
fractional differential equations lead to a well-posed matrix of the resulting 
algebraic system. 



8. Discretization of partial derivatives in time and space 

The simplest implicit discretization scheme for the classical diffusion 
equation is shown in Fig. [2| where the two nodes in time direction are used 
for approximating the time derivative, and the three points in spatial direc- 
tion are used for the symmetric approximation of the the spatial derivative. 
The stencil in Fig. [2] involves therefore only two time layers. If we consider 
fractional-order time derivative, then we have to involve all time levels start- 
ing from the very beginning. This is shown in Fig. |3] for the case of five time 
layers. 

Similarly, if in addition to fractional-order time derivative we also con- 
sider symmetric fractional-order spatial derivatives, then we have to use all 
nodes at the considered time layer from the leftmost to the rightmost spatial 
discretization node. This most general situation is shown in Fig. |4} 

Let us consider the nodes {ih,jT), j = 0,1,2, ...,n, corresponding to 
all time layers at z-th spatial discretization node. It has been shown in 
[51] that all values of a-th order time derivative of u{x,t) at these nodes 
are approximated using the discrete analogue of differentiation of arbitrary 
order: 



{a) (a) 



i,n— 1 



i,2 



%0 



i.n—1 



Ui,2 Ui,l Ui,0 



(31) 
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In order to obtain a simultaneous approximation of a-th order time 
derivative of u{x,t) in all nodes shown in Fig. [T| we need to arrange all 
function values Uij at the discretization nodes to the form of a column vec- 
tor: 



"^nm '^m,n ^m— l,rt • • • ^0,n 

^m,n— 1 ^m— l.n— 1 • • • ^l,ra— 1 ^0,n— 1 



Um,l Um-1,1 ■ ■ ■ %,1 

1 T 

Um,0 Um~l,0 ■ ■ ■ Ml,0 ^0,0 



(32) 



In visual terms of Fig. [T| we first take the nodes of n-th time layer, then 
the nodes of (n — l)-th time layer, and so forth, and put them in this order 
in a vertical column stack. 

The matrix that transforms the vector Unm to the vector uj^"^ of the 
partial fractional derivative of order a with respect to time variable can be 
obtained as a Kronecker product of the matrix which corresponds to 

the fractional ordinary derivative of order a (recall that n is the number of 
time steps), and the unit matrix (recall that m is the number of spatial 
discretization nodes): 

T^) = ® (33) 

This is illustrated in Fig. [5} where the nodes denoted as white and gray are 
used to approximate the fractional-oder time derivative at the node shown 
in gray. 

Similarly, the matrix that transforms the vector U to the vector Ux^^ of 
the partial fractional derivative of order j3 with respect to spatial variable 
can be obtained as a Kronecker product of the unit matrix En (recall that 
n is the number of spatial discretization nodes), and the matrix Rm\ which 
corresponds to a symmetric Riesz ordinary derivative of order f3 [IHl HH] 
(recall that m is the number of time steps): 

S^^l = En® /2(f) (34) 

This is also illustrated in Fig. [5} where the nodes denoted as black and 
gray (corresponding to all discretization nodes from the leftmost to the right- 
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Figure 2: A stencil for 
integer-order derivatives. 



Figure 3: A stencil in case 
of fractional time deriva- 
tive. 



Figure 4: A stencil in case 
of fractional time and spa- 
tial derivatives. 



most one) are used to approximate the symmetric fractional-order Riesz 
derivative at the same node shown in gray. 

Having these approximations for partial fractional derivatives with re- 
spect to both variables, we can immediately discretize the general form of 
the fractional diffusion equation by simply replacing the derivatives with 
their discrete analogs (Fig. |6]). Namely, the equation 

^D>-x^ = /(x,t) (35) 

is discretized as 

®E^-xEn® Ri^^ }u^m = fnm, (36) 

and the matrix of this system has the structure shown in Fig. [7| 

9. Initial and boundary conditions 

Initial and boundary conditions must be equal to zero. If it is not so, 
then an auxiliary unknown function must be introduced, which satisfies the 
zero initial and boundary conditions. In this way, the non-zero initial and 
boundary conditions moves to the right-hand side of the equation for the new 
unknown function. 

10. Implementation in MATLAB 

We provide a set of MATLAB routines for implementing the suggested 
method The function BCRECUR returns the values of the coefficients 
that appear in the fractional difference approximations of fractional deriva- 
tives. The function BAN returns the matrix for the backward difference 
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En ^ Rt 



Figure 5: Discretization of partial 
derivatives. 



d\x\P 



oD^U - D^U = F 



if 



nm J nm 



fn 



Figure 6: Discretization of partial 
derivatives and of the equation 



Figure 7: The structure of the matrix of the resulting algebraic system. 
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approximation of the left-sided ordinary fractional derivative, the function 
FAN returns the matrix for approximating the right-sided ordinary fractional 
derivative, and the functions RANSYM and RANORT return the matrices 



for approximating the symmetric Riesz using the formulas (28) and (29), re 



spectively. The function ELIMINATOR returns the eliminator matrix, and 



the function SHIFT implements the operations (20) and (21). 

The use of these routines is illustrated by the demo functions FRAC- 
DIFFDEMOU, which implements Examples 1 and 2 below, FRACDIFF- 
DEMOY, which implements Examples 3 and 4, and FRACDIFFDEMOY- 
DELAY, which implements Example 5. 

11. Examples 

In this section we introduce several examples illustrating the use of the 
suggested method. 

First, we demonstrate that for the classical integer-order diffusion equa- 
tion our method gives proper results, which are in agreement with the ana- 
lytical and numerical solutions provided in [IB]. 

Second, we obtain the numerical solution of a time-fractional diffusion 
equation. This solution is in perfect agreement with the numerical solution 
obtained in the very recent work [55] by a different approach. 

Then we consider fractional diffusion equation with spatial fractional 
derivative. The fractional derivative with respect to the spatial variable is 
considered as a Riesz fractional derivative. 

After that, we show the results of numerical solution of a general frac- 
tional diffusion equation, where time and spatial derivatives are both of frac- 
tional order - the time fractional derivative is a left-sided Riemann-Liouville 
derivative, and the spatial fractional derivative is a Riesz fractional derivatie. 

Finally, we demonstrate that consideration of partial differential equa- 
tions with fractional derivatives and delays is equally simple in the framework 
of the suggested general approach. 

In all examples, the spatial interval is finite. 

11.1. Example 1: Classical diffusion equation 
Let us start with the classical problem 
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u(0,t) = 0, u(l,t)=0 



(38) 



m(x,0) = 4x(l -x) (39) 

To reduce this problem to a problem with zero initial conditions (the 
boundary conditions are already zero), let us introduce an auxiliary function 



y{x, t) = u{x, t) — u{x, 0) 



(40) 



It follows from (40) and (37)-(39) that the function y{x,t) must satisfy 
dy d'^y 



dt dx"^ 



f{x,t), {with f{x,t) = 8) 



y{0,t) = 0, y{l,t) = 0; 



y{x,0) = 0. 



(41) 



(42) 



The problem (41)-(42) can be discretized using the described method 



see Fig. 6^, which gives 



Er, 



Vn 



fn 



(43) 



where m is the number of spatial discretization intervals and n is the number 
of time steps. 

To obtain the system for finding the unknown values of ynm for the inner 
nodes of the discretization net, we have to use the initial and boundary 
conditions. Since they all are zero, it is sufficient to delete the corresponding 



rows and columns in the system (43), which is easily done with the help of 
eliminators. 

The result of computation of y{x,t) for the spatial step h = 0.1 and the 
time step r = is shown in Fig. M (on the left) for n = 37 time steps. 



These values were chosen for the purpose of comparison with the results from 
[16] . Using (40), we can compute u{x,t), and the result is shown in Fig. [s] 
(on the right). The values of u{x, t) are in perfect agreement with the values 
given in [16] for the same values of h, r, and n. 
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Figure 8: Solutions y{x,t) (left) and u{x,t) (right) of Example 1, with the same values of 
parameters as in [46j. 



11.2. Example 2: Diffusion equation with time fractional derivative 

Now let us consider the equation with the Caputo fractional-order time 
derivative: 



M(0,t) = 0, 



u{l,t) 







(44) 
(45) 



m(x,0) = 4x(l -x) (46) 
Since the Caputo derivative of a constant is zero [2], [50] , for the auxiliary 



function y{x,t) defined by equation (40 ) we obtain a problem with zero initial 



and boundary conditions similar to (41)-(42) 



dx"^ 



f{x,t), {with f{x,t)= 8) 



yiO,t)=0, i/(l,t)=0; 



y{x,0) = 



(47) 



(48) 



This problem can be discretized in the same manner as the previous one 
(refer to Fig. [6]), with the only difference that instead of the first-order time 
derivative we have now a derivative of order a: 



Sl") ® - E„ ® i?!^) 



fn 



(49) 
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where m is the number of spatial discretization intervals and n is the number 
of time steps. 

As above, the use of the zero initial conditions means that the corre- 



sponding rows and columns in the system (49) are removed with the help of 
eliminators. 

The results of computations of y{x, t) and then u{x, t) for a = 1, a = 0.7, 
a = 0.5 with h = 0.05 and r = are shown in Fig. [9| The structure of 
the matrix is the same as shown in Fig. [7| 

Obviously, for a = 1 we have the classical case and the same plots as in 
Fig. |8j and therefore Example 1 is a particular case of Example 2. As a goes 
to zero, the function y{x, t) slowly tends to u{x, 0) = 4x(l — x) for all t. This 
is also not a surprize, because, indeed, for a = the function y{x,t) does 
not depend on t and therefore must satisfy 

y"(x) + 8 = 0, y(0)=y(l) = 0, 

which has the solution y{x) = 4x(l — x). 

It should be noted that almost the same problem as (44)-(46) was nu- 
merically solved in [55] using a very different approach. The initial condition 
in [55] was u{x, 0) = x{l — x). Scaling the plots in figures 1 and 2 in [55] by 
the factor of 4, we obtain the plots which are practically identical with our 
results for u{x,t) shown in Fig. |9} For this comparison we considered the 
shorter interval < t < 0.02 used in 155). 



11.3. Example 3: Dijfusion equation with spatial fractional derivative 

Let us now focus on the role of spatial fractional derivative. For clar- 



ity, let us directly write the following analog of the problem (41)-(42) for 
determining the function y{x,t): 



dy 
dt 



df^y 
d\x\^ 



f{x,t), (with/(x,t) = 8) 



y{0,t) = 0, y{l,t) = 0; 



y{x,0) = 0. 



(50) 



(51) 



where 1 < P < 2. The right-hand side is the same as in (41), but instead 



of second order spatial derivative we deal with the Riesz-Caputo fractional 
derivative. 
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Figure 9: Solutions t) (left column) and u(x, t) (right column) of Example 2, for a = 1 
(top), a = 0.7 (middle) and a = 0.5 (bottom), with spatial step h = 0.05 and time step 
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The problem (50)-(51) can be discretized using the described method 



(52) 



(see Fig. |6|), which gives 



where m is the number of spatial discretization intervals and n is the number 
of time steps, and the corresponding rows and columns in the system (52) 
are removed with the help of eliminators. 

The results of computations for four different values of j3 are shown in 
Fig.flOl 



11.4- Example 4-' General fractional diffusion equation 

Now we can illustrate that the method works also in the case when both 
derivatives are of fractional order. Let us consider the most general situation: 
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Figure 11: Solutions y{x,t) of Example 4, for a — 0.7 and f3 = 1.4 (left), a — 0.7 and 
/3 = 1.8 (middle), a — 0.7 (3 — 2 (right), with spatial step h — 0.05 and time step t — h? /Q. 



QPy 



f{x,t), {with f{x,t)= 8) 



y{0,t) = 0, y{l,t) = 0; 



y{x,0) = 0. 



(53) 



(54) 



The right-hand side is the same as in (41) and (50), but now both deriva- 



tives are allowed to be of non-integer order. 



The problem (53)-(54) can be discretized using the described method 



(see Fig. [6]) , which gives 

®Em- K ® Rl^^}ynn. = fnm (55) 

where m is the number of spatial discretization intervals and n is the number 



of time steps, and the corresponding rows and columns in the system (55) 
are as in all previous examples removed with the help of eliminators. 

The results of computations for some sample combinations of non-integer 



orders a and different values of (3 are shown in Fig. 11 



11.5. Example 5: Fractional diffusion equation with delay 

Finally, let us consider the equation with two Caputo fractional-order 
time derivatives, of which one is with delay 5 (we do not go into the physical 
interpretation of this equation, because physical interpretation of a delayed 
fractional derivative is not known so far, but use it for demonstrating how 
broad can be the field of application of our approach) : 
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\ { ^D^y + ^Dl,y} " |^ = f{x,t) (with f{x,t) ^ 8) (56) 

y{0,t) = 0, y{l,t) = y{x,0) = (57) 

Obviously, for 7 = a and 5 = we have the equation considered in 
Example 2. Let us select the discretization step so that 5 is a multiple of the 



time step t: 6 = kr. Then the problem (56)-(57) can be discretized using 



the described method (see Fig. m^and the equation (21)), which gives 



I ^ (5(-) ® E„ + ® E^) -En® R^£> } yn^ = fnm (58) 



T 



0(7) _ C 771+ r{7) rp+ ql 

+k-Dn — Jn+l,...,n+k J^n+k,k ^n+k ^n+k,k ^1 

where, as above, m is the number of spatial discretization intervals and n is 
the number of time steps, k is the number of time steps corresponding to the 



delay 5, and the appropriate rows and columns in the system (58) are as in 
all previous examples are to removed with the help of eliminators. 

The results of computations for a sample combination of non-integer or- 
ders a, P and 7 and some delays 6 represented by the parameter k are shown 



in Fig. 12 



12. Conclusion 

The suggested method represents a unifying approach to numerical solu- 
tion of partial differential equations of both integer and non-integer order, 
including equations with delays. 

For the sake of clarity, in this article we considered the case of one spatial 
variable. However, the suggested method can be easily extended to the case 
of two and three spatial variables by repeatedly applying the triangular strip 
matrix representations of fractional-order operators in combination with the 
Kronecker matrix product. 

The problems considered in this article are linear, so the resulting systems 
of algebraic equations are linear as well. However, the suggested approach 
can be extended to the case of nonlinear problems, too. 
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Figure 12: Solutions y{x, t) (left column) of Example 5, for a = 0.9, 7 = 0.8, /3 = 1.9, for 
delays 4 = fcr, k = 6, 12, 24, 36. 



The suggested method can be used also for solving partial fractional dif- 
ferential equations appearing from the Laplace equation by replacing second 

order spatial derivatives with fractional Ricsz derivatives. 

The suggested method can be used also for partial fractional FDEs of 
variable and distributed order(s) and for equations with delays. 
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Appendix: sample evaluation of the symmetric Riesz fractional 
derivative 

For (j){x) = x{l — x) and the order of differentiation 1 < (3 < 2 the 
left-sided Riemann-Liouville fractional derivative ^ of the function is 

1-/3 9^2-/3 

Similarly, the right-sided Riemann-Liouville derivative ([T]) of is 



[1-xy-^ 2(1 --^2-/3 

r(2-/3) T( 

Therefore, the symmetric Riesz fractional derivative ([s]) of the function 



[XI is: 



^{o/^f0(a:)+..I)f0(a;)} (61) 

(62) 



a;i-/3 + n-x)^ x^-f^ + (1 - x?-^ 



2T{2-(3) r(3-/5) 
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